[ create a new paste ] login | about

Project: Fast_Fourier_Transform
Link: http://fast_fourier_transform.codepad.org/smayoUBj    [ raw code | output | fork ]

jhafranco - Haskell, pasted on Mar 31:
module Main where
import Data.Complex(Complex((:+)),realPart,imagPart)

main = do
         let z1 = map (roundComplex 5 0.001) $ fft [1,1,1,1,0,0,0,0]
         print z1
         print $ map (roundComplex 5 0.001) $ ifft z1

-- The function 'fft' implements a Fast Fourier Transform (FFT) algorithm for computing
-- the Discrete Fourier Transform (DFT) when N is a power of 2
-- The function 'ifft' implements a Fast Fourier Transform (FFT) algorithm for computing
-- the inverse Discrete Fourier Transform (DFT) when N is a power of 2

fft, ifft :: [Complex Float] -> [Complex Float]
fft = fft' 1
ifft = fft' (-1)

fft' :: Float -> [Complex Float] -> [Complex Float]
fft' e xs = if pow2 (length xs) 
            then let z = sqrt (1 / fromIntegral (length xs))
                  in map ((z :+ 0) *) $ fft'' e xs
            else error "Number of points is not a power of 2"
            where pow2 n
                     | n == 1 || n == 2 = True
                     | otherwise = n `mod` 2 == 0 && pow2 (n `div` 2)

fft'' :: Float -> [Complex Float] -> [Complex Float]
fft'' _ [] = []
fft'' _ [x] = [x]
fft'' e xs = fft'' e (evens xs) <+> t (fft'' e (odds xs))
  where (<+>) r s =  zipWith (+) (r ++ r) (s ++ map negate s)
        evens []  = []
        evens [u] = [u]
        evens (v:_:vs) = v:evens vs
        odds = evens . drop 1
        n = 2 * pi / fromIntegral (length xs)
        t = zipWith (\k z -> z * exp (0 :+ k * n * e)) ([0..] :: [Float])

-- Rounding function
roundComplex :: Int -> Float -> Complex Float -> Complex Float
roundComplex n e x = let zeroFloat e' f = if abs f < e' then 0 else f
                         trim = roundFloat n . zeroFloat e
                         roundFloat m y = fromIntegral (round $ y * 10 ^ m :: Int) / 10 ^^ m
                     in trim (realPart x) :+ trim (imagPart x)


Output:
1
2
[1.41421 :+ 0.0,0.35355 :+ 0.85355,0.0 :+ 0.0,0.35355 :+ 0.14645,0.0 :+ 0.0,0.35355 :+ (-0.14645),0.0 :+ 0.0,0.35355 :+ (-0.85355)]
[0.99999 :+ 0.0,1.0 :+ 0.0,0.99999 :+ 0.0,1.0 :+ 0.0,0.0 :+ 0.0,0.0 :+ 0.0,0.0 :+ 0.0,0.0 :+ 0.0]


Create a new paste based on this one


Comments: